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Abstract 

The objective of this article is to study the behavior of electromagnetic field under 
X-ray diffraction by time-dependent deformed crystals. Derived system of differential 
equations looks like the Takagi equations in the case of non-stationary crystals. This is a 
system of multidimensional first-order hyperbolic equations with complex time-dependent 
coefficients. Efficient difference schemes based on the multicomponent modification of 
the alternating direction method are proposed. The stability and convergence of devised 
schemes are proved. Numerical results are shown for the case of an ideal crystal, a crystal 
heated uniformly according to a linear law and a time-varying bent crystal. Detailed 
numerical studies indicate the importance of consideration even small crystal changes. 
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1. Introduction 

Mathematical modeling of X-ray diffraction by time-dependent deformed crystals refers 
to physical problems of intensive beams passing through crystals. So, a relativistic electron 
beam passes through the crystal target and leads to its heating and deformation. The system 
of differential equations describing X-ray dynamical diffraction by non-stationary deformed 
crystals was obtained in Q. This system looks like the Takagi equations [gj- || in the 
case of non-stationary crystals. Up to now in many ref. (e.g. Q— Q) the theory of X-ray 
dynamical diffraction by stationary crystals for different deformations was developed. Proper 
systems of differential equations are stationary hyperbolic systems for two independent spatial 
variables. In @— @ the solutions of these systems were obtained analytically for some cases 
of deformations. || id devoted to numerical calculation of propagation of X-rays in stationary 
perfect crystals and in a crystal submitted to a thermal gradient. In [||] the theory of time- 
dependent X-ray diffraction by ideal crystals was developed on the basis of Green-function 
formalism for some suppositions. 

The exact analytical solution of the system being studied in this work is difficult if not 
impossible to obtain. That is why we propose difference schemes for numerical solution. To 
solve multidimensional hyperbolic systems it is conventional to use different componentwise 
splitting methods, locally one-dimensional method, alternating direction method and others. 
They all have one common advantage, since they allow to reduce the solving of complicated 
problem to solving of a system of simpler ones. But sometimes they do not give sufficient pre- 
cision of approximation solution under rather wide grid spacings and low solution smoothness 
because the disbalancement of discrete nature causes the violation of discrete analogues of 
conservation laws. The alternating direction method is efficient when solving two-dimensional 
parabolic equations. We use the multicomponent modification of the alternating direction 
method which is devoid of such imperfections. This method provides a complete approx- 
imation. It can be applied for multicomponent decomposition, does not require the operator's 
commut ability. It can be used in solving both stationary and non-stationary problems. 

Difference schemes presented allow the peculiarities of initial system solution behavior. 
The problem of stability and convergence of proposed difference schemes are considered. We 
present results of numerical experiments carried out. We compare efficiency of suggested 
schemes in the case of diffraction by ideal stationary crystal. Tests and results of numerical 
experiments are demonstrated in the case of heated crystal. In our experiments it is assumed 
that the crystal was heated uniformly according to a linear law. The source of crystal heating 
was not specified. It may be the electron beam passing through the crystal. In [Q] we have 
given the formulae which allow to determine the crystal temperature under electron beam 
heating. We present also results of numerical modeling of X-ray diffraction by a time- varying 
bent crystal. The source of crystal bending is not discussed too. 

2. Theoretical Model of X-ray Dynamical Diffraction by 
Time-Dependent Deformed Crystals 

We will use the physical notation jnj. Let a monocrystal plane be affected by some time- 
varying field of forces, which cause the crystal to be deformed. At the same time let a plane 
electromagnetic wave with frequency co and wave vector k be incident on this monocrystal 
plane. We consider two different diffraction geometry which are depicted in Figure. 1. In the 
case of Bragg geometry the diffracted wave leaves the crystal through the same plane that the 
direct wave comes in. In Laue case the diffracted wave leaves the crystal through the back 
plane of the crystal. The electromagnetic field inside the crystal in two-wave approximation 
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a) 



b) 



Figure 1: Diffraction geometry: a) Bragg case, b) Laue case. 



is written in the form: 

D(r, t) = D(r, t) exp(i(fcr - ut)) + D T (r, t) exp(z((fc + r)r - ut)), 

where D and D T are the amplitudes of electromagnetic induction of direct and diffracted 
waves, respectively, and r is the reciprocal lattice vector. 

Let us examine a weakly distorted region in the crystal, where for the deformation vector 
u(r,t) the following inequalities are correct: 



du 



Or 



«1, 



1 du 
c~dt 



« 1, 



where c is the velocity of light. 

We can write :r d = r(l — u). Here T d (r, t) is the reciprocal lattice vector in deformed 
crystal. u(r,t) is the crystal deformation tensor. U{j = l/2{dui/dxj + duj/dxi). Let us call 
considered system of coordinates S. 

To obtain an expansion in series of the reciprocal lattice vector let us pass to a new system 
of coordinates S': r^ = r — u(r,t). Here in each fixed instant of time the Bravais lattice of 
deformed crystal is coincident with one of undistorted crystal in the system S. So, in the 
system S' the crystal structure is periodic. And it is disturbed in the system S. 

Now in S' for electric susceptibility e we can write: 

ef>d;^) = ^2t{i- d ;u;)exp(iT d r d ), 
e(r -u;uj) = Y^ e(r d ; io)exp(ir d (r - it)). 

T d 

Or, finally restoring in the system S, we obtain: 

e(r; u) = e ((! ~ u)t; Lo)exp(ir d r), 
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where e(0; u>) = 1 + g , e((l - u)t;uj) = l+g T (r,t), e(—(l — u)r;u)) = l+g- T (r,t). 

Let us assume that the amplitudes D and D T are changing sufficiently slowly in the space 
and time: 

1 dD 7 



ldD 

k dx; 



< \ D \, 
1 dD 



k dx, 
< \D\, 



< \D T \, 
1 dD 



1, 2, 3. 



< l-Drl- 



w <9i ' uj dt T ' 

Then from Maxwell's equations the following system of differential equations was derived 



2idD 2i, Jn n n 

+ -3 fcgrad£> + xoD + X rD T = 0, 



uJ 



dt k 2 



2i dD 2i 

tt: 1 + j^krgvadD-r + X-r D + (xo ~ ol(t, t) - s(r, t))D 7 

uj at k z 



(2.1) 



where 



, , 2fc T grad(ru) (r 2 + 2kr) 
a(r, t) = a p , a = p 



. , 2 / du\ 



c 



Xo 5 X±t are the zero and ±r Fourier components of the crystal electric susceptibility. 

The difference between our system (2.1) and the Takagi equations ||] is in the term a, 
which depends on time now, and in the appearance of the term s. 

Let us rewrite the system (2.1) in the generalized form having picked out vectors of a- 
polarization from amplitudes of electromagnetic induction D and D T and having specified 
three independent variables t, z, x. The spatial variable y is a parameter. One can write a 
full three-dimensional system. 

dD A dD A dD n „ n n 
— + A n — + A u — + Q n D + Q 12 D T = 0, 
at dz dx 



dD T dD T dD T 

+ A 21 ^+ A 22 ^ + Q21D + Q 22 D T = 0, 



dt 



dz 



dx 



where 



ck~ 



111 



A 12 = 

k 



A 



21 



A 22 = CK '- 



Q 



11 



k ' k 

0.5iu)Xo, Q12 = -0.5iwxr, 



(2.2) 



(2.3) 



Q21 = -OMujx-t, Q22 = -0.5iuj(xo - ot(z, x, t) - s(z, x, t)). 



(2.4) 



Initial and boundary conditions are written in the domain G = {(z, x, t), < z < L z , < 
x < L x ,0 < t < T}. In the Bragg case the boundary conditions are written as follows: 



D(Q,x,t) = D 0l 
D T (L z ,x,t) = 0, < x < L x , t>0. 



(2.5) 



In Laue geometry, where the diffracted wave leaves the crystal through the crystal back 
plane, the boundary condition for amplitude D r should be written at z = 0. 
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As is known [ 11 1 , the exact solution of the stationary X-ray diffraction problem has the 
following form: 

D = c\ exp(ik5iz) + C2 exp(ik52z), 

D T = c\S\ exp(ifcJiz) + C2S2 exp(ik52z), (2-6) 
where 61,62 are the solutions of the dispersion equation: 

(2^7o - Xo)(257i - «o - Xo) - XrX-r = 0; 
_ 2^70 - xo ■ , 9 . 

Xr 

70 and 71 are the cosines of the angles between k and k T , respectively, and the z axis; 

-D s 2 e2(L z ) Z? siei(L 2 ) 
c i = rr^ rrr> c 2 ~ 



siex(L z ) - s 2 e 2 (L z ) siei(L z ) - s 2 e 2 (L z ) 

Here and below the following designations are used: 

e\{z) = exp(ik5iz), e2(z) = exp(ik62z). 

Let us impose initial conditions corresponding to the exact solution of the stationary 
X-ray diffraction problem in an ideal crystal: 

D(z,x,0) = c\ei(z) + C2e 2 (z), 

D T (z,x,0) = c\Sie\{z) + 0252^2(2), < z < L z , < x < L x . 

In the X-ray range the amplitudes (2.6) oscillate with sufficiently high frequency. For 
large thickness of crystal it is complicated to obtain good numerical solutions of system (2.2) 
with coefficients (2.4). So, let us find solution of (2.2) for functions D(z,x,t) and D T (z,x,t) 
which vary more slowly than e\(z) or e2{z): 

D(z,x,t) = D(z,x,t)(aei(z) + C2e2(z)), 
D T (z,x,t) = B T (z,x,t)(ciSiei(z) + c 2 s 2 e 2 (2;)). (2.7) 

Then the coefficients (2.4) have to be presented by the formulae: 

Xo - 1k z jk{c\6\e\(z) + C25 2 e2(z)) 



Qu = — 0.5iu; : 

Q12 = -0.5iojx 
Q21 = -0.5iux- 



ciei(z) + c 2 e 2 (2) 

cisiei(z) + c 2 s 2 e2(^) 
ciei(z) + C2e 2 (z) 

c 1 e 1 (z) + c 2 e 2 (z) 
r cisiei(z) + c 2 S2e 2 (z) ' 



^ nr;- Xo - oi{z, x, t) - s(z, x, t) - 2k TZ /k(cisi6 1 e 1 (z) + C2S 2 S 2 e2(z)) 

V22 — —u.btu -— — . (/.8) 

cisiei(z) + C2S2e 2 {z) 

The boundary conditions (2.5) take the form: 

D(0,x,t) = 1, 

D T (L z ,x,t) = 1, 0<x<L x , t>0. (2.9) 
For this case the initial conditions have to be equal to 1 too. 



5 



3. Numerical Analysis 

The original differential problem is a system of multidimensional first-order differential equa- 
tions of hyperbolic type with complex-valued time-dependent coefficients. The numerical 
schemes employed in this work are based on the multicomponent modification of the alter- 



nating direction method. This method was originally developed in [10|. It turned out to be 
an effective way for efficient implementations of difference schemes. This method is econom- 
ical and unconditionally stable without stabilizing corrections for any dimension problems of 
mathematical physics. It does not require spatial operator's commutability for the validity 
of the stability conditions. This method is efficient for operation with complex arithmetic. 
Its main idea is in the reduction of the initial problem to consecutive or parallel solution of 
weakly held subproblems in subregions with simpler structure. That is why it allows us to 
perform computations on parallel computers. The main feature of this method implies that 
the grids for different directions can be chosen independently and for different components 
of approximate solution one can use proper methods. 

Let introduce the Hilbert space of complex vector functions: H = L 2 (G). In this space 
the inner product and the norm are applied in the usual way: 

(u,v) = / u(x)v(x)dx, \\u\\ = (u, u) 1 / 2 . 



G 

In H our system (2.2) is hyperbolic. 

3.1. Efficient Schemes for Solving Hyperbolic System in Two Space Dimensions 

We use the following notation |[12| : 

Ux = (Ui+i — Ui)/h x — right difference derivative, 
m = (Vi ~ Vi-i)/h x — left one, y-i = y(xi); 

yt = (y-y)/ht, y = y(tk+i), y = y(tk)- 

Let us replace the domain G of the continuous change of variables by the grid domain 
G zxt = {(zi,Xj,t k );zi = ih z ,i = 0, 1, . . . , iVi, iVi = [L z /h z ],Xj = jh x ,j = 0, 1, . . . , N 2 , 

N 2 = [L x /h x ],t k = kht, k = 0, 1, . . . , N 3 , N 3 = [T/h t ]}. 

The following system of difference equations approximates on G zx t the system (2.2) with 
coefficients (2.3)-(2.4): 

D] + A n Dl + A l2 Dl + QnD 1 * + Qi 2 D l T * = 0, 
D\ t + A 21 D l TZ + A 22 D 2 TX + Q 21 D l * + (Q^)* = 0, (3.10) 



D 2 t + A XX B\ + A 12 D 2 X + QnD 1 * + Q l2 D l T * = 0, 
D 2 Tt + A 2l D l TZ + A 22 D 2 TX + Q^D 1 * + (Q^DD* = 0, (3.11) 

where D* = 0.5(A + A-i), Q22* = 0.5(Q 2 2(^i-l, Xj, tk+i) + Q22(Zi,Xj,t k+ i)) for the first 
equations of (3.10) and (3.11) and D* = 0.5(Di + -D«+i), Q 22 = 0.5(Q 22 (zi, Xj, t k +i) + 
Q22(zi+i,Xj,tk+i)) for the last ones. 
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For coefficients (2.3), (2.8) the system (2.2) can be approximated by the system of differ- 
ence equations of the following form: 

D\ + A XX D\ + A 12 D 2 X + QnD 1 + Q l2 D l T = 0, 
Dlt + A 21 D l TZ + A 22 D 2 TX + Q 2l D l + Q^D\ = 0; 

d 2 + And* + A 12 b 2 x + Qub 1 + Q X2 b\ = o, 
D 2 t + A 2X b\ z + ^22^ + Q 2 ib l + q^dJ: = o. (3.i2) 

In cited schemes the directions of difference derivatives with respect to x (left or right) 
are selected in dependence of waves directions. D , D 2 , D\ and D 2 are two components 
of approximate solutions for D and D T , respectively. One can choose any of these two 
components or its half-sum as a solution of (2.2). The boundary and initial conditions are 
approximated in the accurate form. 

In Laue case where the diffracted wave moves on a positive direction of the z axis, for D T 
we should write left difference derivatives with respect to z. 

The schemes (3.10)— (3.11) and (3.12) are completely consistent. The consistency clearly 
follows from the manner in which these schemes were constructed. On sufficiently smooth 
solutions they are of the first order approximation with respect to time and space. We can 
give the difference scheme of the second order approximation with respect to z. In this case 
it should be rewritten (3.10): 

D]* + A 1X D\ + A 12 D 2 X + + Q 12 Dl* = 0, 

Dl* + A 2l D l TZ + A 22 D 2 TX + Q 2l D l * + (Q^DD* = 0. (3.13) 

The scheme for the second component (3.11) is not changed. 

One can write a scheme of the second order approximation with respect to time. But 
as has been shown in numerical experiments it does not lead to sensible changes in solution 
pattern. 

For the difference schemes presented, the stability relative to initial data and also the 
convergence of the difference problem solution to the solution of differential problem (2.2) 
can be proved. This follows from the properties of the multicomponent modification of the 



alternating direction method [10|. Let us prove the corresponding Theorems. 



3.2. Stability and Convergence of Difference Schemes 

We use the energy inequalities method |0|. Let rewrite the system (3.10)— (3.11) in the 
form: 

D f 1 + A 1 (D 1 )+A 2 (D 2 ) =0, (3.14) 
D t 2 + A 1 (D 1 )+A 2 (D 2 ) = 0, (3.15) 

where 

□=(£), A 1(D) = ( A M t:^:^uY ) • « D > = ( A 2 D 

Let us introduce the following notation: 

y' = Re{y), y" = Im(y), y = y'-iy". 
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We use the inner products: 

N-l 

(y,v)u = E h Vi v ^ 
i=i 

where u = {xi = ih,i = 0, 1, N, Nh = L} is a one-dimensional grid; 

Ni-l N 2 -l 

(Y,V) Gzxt = (Y,V)= E 

i=i j=i 

In addition to this let us introduce the norm: 

\\y\\ = ^y), |y| 2 = ||y'|| 2 + ||y"|| 2 . 

Lemma. // y(x ) = then (y,yx)u> > 0, if y{x N ) = then (y,y x )u> < 0. 

Proof. Let us write the following transformation chain: 

N-l N-l 

{y,yx)u = E ViiVi ~ fi-l) = E (fi ~ fi-l ~ fi-l(f< ~ fi-l)) = 
i=l i=l 
JV-1 JV-2 

yjv_l - 2/o - E ~ W-i) = VN-l - 2/0 - E ~ Pi) = 

i=l i=0 
N-2 N-2 N-2 N-2 

y 2 N-i - yl + E y 2 - °- 5 E (w+i + y*) 2 + °- 5 E y 2 +i + °- 5 E yt = 

i=0 i=0 i=0 i=0 

AT-2 JV-2 N-2 

0.5y 2 N -i ~ 0-5y 2 + E ^ + E S/i+i ~ °- 5 E + ^) 2 - 
i=0 i=0 i=0 

The first term in the last expression is greater than or equal to 0, the second is equal to 
0. Taking into consideration the inequality: 

(yi+i+Vi) 2 <2(y 2 +1 +y 2 ), 

we obtain: (y, yx)u> > 0. Second Lemma's inequality is proved similarly. □ 

Teorem 3.1. The difference scheme (S.lO)-(S.ll) is unconditionally stable relative to 
the initial data. For its solution the following estimates hold: 

|Df < M (|D»(*o)| 2 + IMD 1 ^)) + A 2 (D 2 (t ))| 2 ) , (3.16) 

where M is a bounded positive constant independent of grid spacings, i = 1,2. 
Proof. Multiply (3.14) by (A^D 1 ))^ (3.15) by (A 2 (D 2 )) t and sum : 

((D t 1 )',(A 1 (D 1 ));) + ((D t 1 r,(A 1 (D 1 ));') + 

((D 2 )',(A 2 (D 2 ));) + ((D 2 )",(A 2 (D 2 ));') + 

((Ai(D 1 )) / , (Ai(D 1 ))^ + ((Ai(D 1 ))", (Ai(D 1 ))j'j + 

((A 2 (D 2 ))',(A 1 (D 1 ));) + ((A 2 (D 2 ))",(A 1 (D 1 ));') + 

((A 1 (D 1 ))',(A 2 (D 2 )X) + ((A 1 (D 1 ))",(A 2 (D 2 ));') + 

((A 2 (D 2 ))', (A 2 (D 2 ))^) + ((A 2 (D 2 ))", (A 2 (D 2 ))?) = 0. (3.17) 
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Let us multiply (3.17) by h t and take into account the form of time derivatives. Then we 
obtain: 

$(D) + 0.5 (jAitD 1 ) + A 2 (D 2 )| 2 + lA^D 1 ) + A 2 (D 2 )| 2 ) + 
0.5 (IMD 1 )! 2 + |A 2 (D 2 )| 2 + lA^D 1 )! 2 + |A 2 (D 2 )| 2 ) - 
((Ai(D 1 )) / , (Ai(D 1 ))'^ - ((A 2 (D 2 ))',(A 2 (D 2 ))') - 
((MD 1 ))", (A^D 1 ))") - ((A 2 (D 2 ))", (A 2 (D 2 ))") = 0, (3.18) 

where 

*(D) = ht ((DlY, (Ai(D x ))^ + h t ((A 1 )", (Ai(D 1 ))^ + 

h t ((A 2 )', (A 2 (D 2 ))^) + ht ((A 2 )", (A 2 (D 2 ))?) . (3.19) 

Let us re-arrange (3.18): 

$(D) + O.SIA^D 1 ) + A 2 (D 2 )| 2 - O^A^D 1 ) + A 2 (D 2 )| 2 + 
0.5/if |(Ai(D 1 )) t | 2 + 0.5/i 2 |(A 2 (D 2 )) t | 2 = 0. 

Or, more properly, 

$(D) + O.SIA^D 1 ) + A 2 (D 2 )| 2 < O^jA^D 1 ) + A 2 (D 2 )| 2 . (3.20) 

If we prove that ^(D) > then the estimate (3.16) will be obtained for i = 2. We consider 
the first term in (3.19) 

((A 1 ) / ,(Ai(D 1 )l) + ((A 1 r,(Ai(D 1 )?) = 
((Dl)',(A 11 Dl)' t ) - ((£$)' \Q 11 (D 1 *)?) - ((Dj)',Q 12 (D l T *y;) + 
((£&)', (A 21 D 1 TZ )' t ) - ((D 1 Tt )',Q 21 (D 1 *)!) - ([Di t y,{Q 22 Di*yi) + 
((Dl)",(A 11 Dl)l) + ((Dl)",Q 11 (D 1 *y t ) + {(Dl)",Q l2 (Dl*y) + 
{(Dl t )",(A 2l D\ z yi) + {{D\ t )\Q 21 {D^y t ) + ((Dl t )",(Q 22 Dl*y t ) > 0, (3.21) 

as since from Lemma we have: 

((A 1 )'. (^n^)O > 0, ((A 1 )". (AuD^) > 0, A 11 = (A 11 )' > 0; 

[{D\ t y, {A 21 D\ z y t ) > o, ((Di t y, (A 2l D l TZ yi) > o, a 21 = (A 21 y < o. 

In (3.21) we took into account that the coefficients Q from (2.4) were pure imaginary, Q\ 2 = 
Q 2 \. Also we have written out corresponding inner products. 

In the same fashion as Lemma we obtain for the second term in (3.19): 

((A 2 )', (A 2 (D 2 X) + ((A 2 )", (A 2 (D 2 );') > 0. (3.22) 

So, expressions (3.21) and (3.22) mean that operators Ai and A 2 are positive definite. 
If we repeatedly apply the expression (3.20) to the right and if we take into account that 

D 2 = -A 1 (D 1 )-A 2 (D 2 ), (3.23) 
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we have the estimate (3.16) for i = 2. 

Now let us obtain the estimate (3.16) for i = 1. Multiply D 1 by D 1 and take into account 
the expression which follows from our designations: 



D^D 1 -^ (A 1 (D 1 ) + A 2 (D 2 ) 
Using e-inequality fl2| |(it,v)| < e||u|| 2 + l/(4e)||w|| 2 , (e > 0), we have : 



\2\l2 



ID 1 ! 2 = [D 1 - ht (AiiD 1 ) + A 2 (D 2 )) , Di - ht (Ai(Di) + A 2 (D2) 

ID 1 ] 2 + ^lAiCD 1 ) + A 2 (D 2 )| 2 - 2h t ((D 1 )', (^(D 1 ) + A 2 (D 2 ))'^ 
-2h t ((D 1 )", (A^D 1 ) + A 2 (D 2 ))") < MxlD 1 ! 2 + M 2 |A 1 (D 1 ) + A 2 (D 2 )| 2 . 

Let us consider the second term in previous inequality: 

IMD 1 ) + A 2 (D 2 )| 2 = IMD 1 ) + A 2 (D 2 ) - h t (A 2 (D 2 )) t \ 2 
IMD 1 ) + A 2 (D 2 )| 2 + / l 2 |(A 2 (D 2 )) 4 | 2 - 2ht ((AiiD 1 ) + A 2 (D 2 ))', (A 2 (D 2 ))^ 

2ht ((A 1 (D 1 )+A 2 (D 2 ))",(A 2 (D 2 ))^) < |A i (D 1 ) + A 2 (D 2 
This was obtained by taking into account (3.23) and (3.22). Finally we have: 

ID 1 ] 2 < M 3 (ID 1 ! 2 + IAiCD 1 ) + A 2 (D 2 )[ 2 ) . 

Repeated application of this inequality to the right yields the estimate (3.16) for i = 1. □ 

We denote the discretization error as Z l = D* — D, i = 1, 2, where D is the exact solution 
of initial differential problem. 

Teorem 3.2. Let the differential problem (2.2)-(2.5) have a unique solution. Then 
the solution of the difference problem (3.10)-(3.11) converges to the solution of the initial 
differential problem as ht, h z ,h x — > 0. The discretization error may be written as 

\Z l \ < 0(h t + h z + h x ). 

Proof follows immediately from consistency of the scheme (3.10)-(3.11), Theorem 3.1 
and Lax's Equivalence Theorem [13]. □ 

The stability and convergence of schemes (3.12), (3.13)-(3.11) can be proved in an anal- 
ogous way. 

4. Results of Numerical Experiments 

4.1. Diffraction by Ideal Crystal 

The problem of studying of electromagnetic fields under X-ray diffraction inside the crystal 
target is constituent of the problem of modeling intensive beams passing through crystals, X- 
ray free electron laser and others. Therefore let analyze the operation of schemes (3.10)-(3.11) 
and (3.12) in the case of ideal absorbing crystals. 

Figure 2-5 display results of numerical experiments in the crystal of LiH. This crystal 
was chosen because of small absorption. The design parameters were the following. The 
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Figure 2: Amplitudes of direct wave for L = 0.1 cm, a) N\ = 100, b) N\ 



200. 



frequency u> was equal to 1.7 • 10 19 sec _1 . The diffraction plane indexes were (2, 2, 0). The 
Bragg angle was equal to 0.39. The angle between the direct wave vector and the z axis was 
equal to 0.83. This case corresponds to the total internal reflection region. We compare the 
numerical results obtained with the analytical solution (2.6). On our figures this solution is 
depicted by red curves. So, Figure 2 presents numerical results of the scheme (3.10)-(3.11) 
for N\ = 100 and Ni = 200 respectively. From these plots we notice that the grid dimension 
N\ = 100 gives good agreement with analytical results and N\ = 200 gives ideal agreement. 
But for greater thickness of L z = 0.3 cm only N\ = 200 gives a more or less acceptable fit, 
as is obvious from Figure 3. 

Thus, for large thickness of crystal we will use the scheme (3.12) with coefficients (2.8). 
It is evident that the amplitudes D and D r from (2.7) should be equal to 1 after installation 
of a stationary regime in the system in the case of an ideal crystal. Figure 4 demonstrates 
the distinction between numerical results for N\ = 20 (curve 1), N\ = 50 (curve 2) and 1. 
When Ni = 100 the agreement is ideal. For the amplitudes of diffracted waves we can show 
similar figures. Figure 4 b) presents numerical and analytical solutions for the amplitudes of 
direct wave when Aq = 50. 

Let us clear up how the scheme (3.12) functions with coefficients (2.3)-(2.4). Figure 5 
gives an idea that in the case of oscillations of amplitudes this scheme operates badly even 
for small thickness of crystal. 

Let us show numerical results for the crystal of Si with small thickness L z = 0.005 
cm. They are more visual because the absorption coefficient of Si is large. We have used the 
following geometry parameters: the diffraction plane indexes are (2, 2, 0), uj = 6.9 • 10 18 sec _1 . 
The Bragg diffraction case was modeled with the Bragg angle equal to 7r/4. Figure 6 depicts 
curves of amplitudes of direct and diffracted waves in comparison with analytical solution in 
the ideal crystal (2.6) (red curves). Figure 7, a) represents comparison between results of 
schemes (3.10)-(3.11) and (3.12) with coefficients (2.4). As stated above, the simplest scheme 
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Figure 3: Amplitudes of direct wave for L = 0.3 cm, a) N\ = 100, b) N\ = 200. 
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Figure 4: a) Reduced amplitudes of direct wave for L = 0.3 cm, 1 — iVi = 20, 2 - 
A^i = 50 by scheme (3.12) with coefficients (2.8). b) Amplitudes of direct wave for L = 0.3 
cm, N\ = 50 (numerical and analytical solutions). 
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Figure 5: Amplitudes of direct wave for a) L = 0.01 cm, 1 — Ni = 100, 2 — Aq = 200 b) 
L = 0.02 cm, 1 — Ni = 100, 2 — Aq = 200 by scheme (3.12) with coefficients (2.4). 

does not work well in our conditions. Figure 7, b) shows the behavior of two components 
D 1 and D 2 of numerical solution. As may be seen from this figure, both of two components 
converge well to the analytical solution. That supports once again our statements given 
above. 

4.2. Diffraction by Time-Dependent Heated Crystal 

Let us compare obtained numerical results for the crystal heated to the temperature 
T = T — Tq K with analytical solutions of the stationary linear X-ray diffraction problem in 
crystal heated to T K. To is an initial temperature of the crystal. The source of crystal heating 
and deformation was not specified. It was supposed that the crystal was heated uniformly 
according to the linear law: T(t) = Tq + at, where a is the rate of heating. The data for 
the stationary linear X-ray diffraction problem (xo-i X±r(T)) were obtained from the program 
flijl . Let emphasize that only the values xo and X±t{Tq) in numerical calculations by (3.10)- 
(3.11) are needed. They can be taken from reference books or from the program JH| ]. While 
computing the values of a and s from the coefficient Q22 are recalculated depending on the 
variation of the deformation vector u(z,x,t). 

We demonstrate numerical results for the crystal of Si. Figure 8 depicts curves of ampli- 
tudes of diffracted wave in the crystal of Si for the above parameters and T = 10 K (curves 
1), T = 15 K (curves 2), T = 20 K (curves 3), respectively. The initial temperature To was 
equal to 293 K. The heating rate a was equal to 5 • 10 10 K/sec. The curves of each of pairs 
of curves correspond to the numerical solution of the X-ray dynamical diffraction problem 
in time- varying heated crystal and to the analytical solution of the stationary linear X-ray 
diffraction problem in the heated deformed crystal. Figure 9 and Figure 10 show evolution of 
direct and diffracted wave amplitudes under crystal heating to 100 K, respectively. Appar- 
ently, up to T = 16 K the modulus of the amplitude of the diffracted wave coming out the 
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Figure 6: Amplitudes of a) direct wave and b) diffracted wave in the crystal of Si by scheme 
(3.10)-(3.11) with coefficients (2.4). 




Figure 7: Amplitudes of direct wave a) for the schemes (3.10)-(3.11) (curve 1) and (3.12) 
(curve 2); b) two components D 1 (curve 1) and D 2 (curve 2) of numerical solution (3.10)- 
(3.11). 
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Figure 8: Amplitudes of diffracted wave when crystal of Si heating to T = 10 K (curves 1), 
T = 15 K (curves 2), T = 20 K (curves 3). 





Figure 10: Evolution of diffracted wave amplitude under crystal heating to 100 K. 
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Figure 11: Model of crystal bending. 

crystal at z = decreases abruptly. This can be explain by the fact that under heating the 
parameter a of deviation from the exact Bragg condition increases and diffraction disrupts. 
Such an analysis of results demonstrates that proposed mathematical model and effective 
numerical algorithm allow to obtain distributions of electromagnetic waves amplitudes in 
non-stationary crystals with sufficient precision. 

4.3. Diffraction by Time-dependent Bent Crystal 

Let us examine the following model of bent in time crystal. As before, we do not specify 
the nature of bending (mechanic, temperature or other). We suppose that the crystal is bent 
according to the law 

u z (x,t) = atx 2 (4.24) 

(see Figure 11), where a is the rate of bending. (4.24) was obtained from the following 
assumptions. The crystal plane formula at the point zq is: 

Z = Zq. (4.25) 

The parabola formula at the point zq has the form: 

atx 2 = z - z . (4.26) 

The difference between (4.25) and (4.24) is the z component of the deformation vector 
u(z,x,t). The component u x (z,x,t) can be derived from the formula: 

u x (z,x,t) = x - x, 

where xq is found from the curve distance formula 

J sfl + 4a 2 t 2 £ 2 d£ = x, 
o 
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Figure 12: Amplitudes of diffracted wave under crystal bending at x = 0.1 cm. 



or: 

y ^1 + Aa?t 2 xl + — ln| 2atx + + 4a 2 t 2 sg 

But the estimations show that when the magnitude of u z is not large (of the order 10 -6 ), the 
magnitude of u x is of the order 10~ 8 and can be neglected. 

We have realized numerical experiments to find out how crystal bending affects the diffrac- 
tion pattern. We consider the crystal of Si with the rate of bending a = —2.5 • 10 6 (cmsec) _1 . 
Figure 12 illustrates evolution of diffraction at the point x = 0.1 cm. The pattern of distri- 
bution of electromagnetic field is similar to one under crystal heating. When moving to the 
central point x = of the crystal the magnitude of bending becomes smaller. The diffraction 
pattern should be less changed. This fact was confirmed during numerical experiments. 

We considered simple models of crystal heating and bending. More accurate ones can 
be taken, for example, from [l5| . In Q-0 a general dynamical theory of X-ray diffraction 
from a homogeneously bent stationary crystals was developed. But their analytical formulae 
are complicated enough. That is why the analysis of numerical results obtained from our 
program and their analytical results will be the aim of another paper. 

5. Summary 

Presented difference schemes and numerical algorithms allow to examine waves amplitudes 
evolution in non-stationary crystals with sufficiently precision. Numerical calculations show 
that even small non-stationary crystal deformations lead to considerable changes in the 
diffraction pattern. So, mathematical model and numerical method presented can be used in 
mathematical modeling of intensive beams passing through crystals. 
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